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ABSTRACT 

Recent observations show that the coohng flows in the central regions of galaxy clusters 
are highly suppressed. Observed AGN-induced cavities/bubbles are a leading candi- 
date for suppressing cooling, usually via some form of mechanical heating. At the 
same time, observed X-ray cavities and synchrotron emission point toward a signifi- 
cant non-thermal particle population. Previous studies have focused on the dynamical 
effects of cosmic-ray pressure support, but none have built successful models in which 
cosmic-ray heating is significant. Here we investigate a new model of AGN heating, in 
which the intracluster medium is efficiently heated by cosmic-rays, which are injected 
into the ICM through diffusion or the shredding of the bubbles by Rayleigh- Taylor 
or Kelvin-Helmholtz instabilities. We include thermal conduction as well. Using nu- 
merical simulations, we show that the cooling catastrophe is efficiently suppressed. 
The cluster quickly relaxes to a quasi-equilibrium state with a highly reduced accre- 
tion rate and temperature and density profiles which match observations. Unlike the 
conduction-only case, no fine-tuning of the Spitzer conduction suppression factor / 
is needed. The cosmic ray pressure, Pc/Pg ^ 0.1 and VPc ^ O.lpg, is well within 
observational bounds. Cosmic ray heating is a very attractive alternative to mechan- 
ical heating, and may become particularly compelling if GLAST detects the 7-ray 
signature of cosmic-rays in clusters. 

Key words: cooling flows - galaxies: clusters: general - cosmic rays - instabilities - 
X-rays: galaxies: clusters 



1 INTRODUCTION 

Clusters of galaxies contain a large amount of hot diffuse gas 
which emits prolifically in thermal X-rays. The X-ray surface 
brightness of many galaxy clusters is strongly peaked in the 
central regions, where the cooling time is much shorter than 
the Hubble time. In the absence of any heating sources, the 
radiative cooling due to this emission will induce a subsonic 
inflow of gas in order to maintain pressure equilibrium, lead- 
ing to substantial dropout of cold gas in the inner regions 
of rich clusters. The mass deposition rates at the cooling ra- 
dius, where the cooling time equals to the age of the cluster, 
were estimated to be as much as several hundred Mq yr~^ 
in some clusters (see [Pabian 199^, for a review). However, 
recent high-resolution Chandra and XMM-Newton observa- 
tions show that while the temperature is declining toward 
the center, there is a remarkable lack of emission lines from 
the gas at temperature belo w about ~ 1/3 of the amb i- 
ent cluster temperat ure (e.g., [Peterson et al.l (1200 ll. |2003|): 
iTamura et al] (I2OOII ): for a review see iPeterson fc FabiarJ 
1 20061 )'). Moreover, the spectroscopically determined mass 



deposition rates are about 10 times smaller than the clas- 
sic values es timated from X-ray lu minosity within the cool- 
ing regions l|Voigt fc Fabianll2004 ). and thus limited to at 
most ~ few X lOM0yr~^. These discrepancies suggest that 
mass dropout is prevented, or significantly reduced, by heat- 
ing sources. Alternatively, the gas cools with line emis- 
sion highly suppressed by mixing, inhomogeneous metallic- 
ity distributions, differen t ial absorption, or photoionizatio n 
(|Fabian et al.ll200ll . |2002| : iMorris fc Fabiaiill2003l : Iohll2004l ) , 
though a fully successful model which explains the observa- 
tions on such grounds is still lacking. 

Many mechanisms for heating the intracluster gas have 
been put forth recently, including transport of heat from the 
hot outer regions of the cluster by thermal conduction (e.g. 
tZakamska & Narayan 2003: Voigt fc F abian 200^ or tur- 
bulen t mixing (|Kim fc' Naravaiill2003bi : lDennis fc ChandrarJ 
I2OO5I ) and heating by outflows, bubbl es or sound waves from 
a central active ga l actic n uclei (e.g., Churazov et al.l 



Kaiser 200d; iRuszkowski fc BegelmaiJ 

2004). Recent theoretical and numerical 
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work (e.g. [Naravan fc Medvedevll200ll : ICho et al.|[20o3 ) has 
shown that a turbulent magnetic fleld may not be as efficient 
in suppres sing thermal conduction as pr eviously thought. In 
particular. iNarayan fc MedvedevI (|200ll ) showed that the ef- 
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fective thermal conductivity k in a turbulent MHD medium 
is a substantial fraction (~ 1/5) of the classical Spitzer value 
KSp if magnetic turbulence extends over at least two decades 
in scale. However. Voigt fc Fabian ( 2004.) found that thermal 
conduction alone is usually insufficient to heat the gas in the 
inner parts of hot clusters and most regions of cool clusters. 
Moreover, balancing cooling with t hermal conduction alone 
is generally globally unstable (e.g. iKim fc NaravanI l2003al ) 
and, as we shall see in § 14.31 requires fine-tuning. 

On the other hand, heating from a central AGN pro- 
vides a self-regulating feedback mechanism, which may play 
a key role in halting a cooling catastrophe in the intra-cluster 
medium (ICM). Here the cooling flow triggers AGN activ- 
ity and heating, which in turn suppresses the cooling catas- 
trophe. The flow thus automatically adjusts itself to a low 
value of the accretion rate, depending mainly on the feed- 
back coefficient (e in eq. [H]), so there is no fine-tuning 
problem. Recent observations suggest that active galactic 
nuclei (AGNs) at the centers of the clusters may interact 
with and substantially heat the ICM. Radio-loud a ctivity 
is ver y common at the centers of cool core clusters (|Burnsl 
Il99d ). Recent high-resolution X-ray observations have also 
revealed cavities or bubbles h aving sizes of a few kiloparsecs 
in many galaxy clusters (e.g. iBirzan et al.ll2004h . 

Such AGN feedback heating models have received a 
good deal of attention in recent years. Numerical sim- 
ulations of AGN bubbles (and jets, in some cases) in 
the ICM have been performed by a number of au- 
thors (e.g . 'BriiEEcn & Kaisei"2002'; 'Ruszkowski et al.l l2004 
IVernaleo & Rcvnolds 2006; Revnolds et al. 2005). These 
simulations usually focus on the heating of the ICM by the 
pdV work of the expanding bubbles, viscous dissipation of 
emitted sound waves or mixing of the hot bubble plasma 
with the ICM. Most simulations assume that these X-ray 
cavities are filled with low-density gas at very high temper- 
ature (e.g., ~ 100 keV). However, radio synchrotron and in- 
verse Compton emission has been observed from many cav- 
ities, suggesting a significant non-thermal component, such 
as cosmic rays and magnetic fields, in these cavities. Cos- 
mic rays (CRs) can indeed be injected at the tip of a ra- 
dio jet, which moves supersonically in the ICM at its initial 
momentum-driven phase and forms radio cocoons or bub- 
bles at its later stage. A substantial amount of cosmic rays 
may then es cape from these buoyantly rising bubbles (e.g. 
EnfiUnI I2OO3I) and heat the ICM. X-ray observations (e.g. 



Fabian et al.ll2006l ) show that some bubbles do remain sta- 



ble even far from cluster centers, but a significant fraction 
of the bubble population could be shredded or disrupted as 
they rise through the ICM. Indeed, numerical simulations 
show that the surfaces of buoyant bubbles are highly suscep- 
tible to disruptio n by Rayleigh- Taylor and Kelvin-Helmholz 
instabilities fe.g. iBriiggen fc Kaiseij [200^ ) . unless viscosity 
or magnetic fields are invoked to suppress these in stabilities 
fe.g. lRevnolds et al.ll2005l : lRuszkowski et al.ll2007h . Hence, it 
is of great interest to study the consequences — particularly 
the heating effect on the ICM — of cosmic rays leaked or dis- 
rupted from radio cocoons or bubbles. 

Cosmic rays may also be produced by other processes 
near the central AGN. Structure formation shocks, merger 
shocks and s upernovae may also inject cosmic rays into 
the ICM (e.g. IVolk et~al]|l996l : [Berezinskv et al.lll997l ). Di- 
rect evidence for the presence of an extensive population of 



nonthermal particles in galaxy clusters comes from the ob- 
servation of diffuse radio synchrotron emission (e.g., radio 
haloes, mini-haloes a n d relics) in many massive clusters (e.g. 
Brunetti et"all I2OOII : IPetrosianI I2OOII : IPfrommer fc Eni3Unl 
20041 : lFerettill2005l ): recent Chandra and XMM observations 



also show evidence for a s ignificant nonthermal particle pop- 
ulation within the ICM JSanders et all l2005l : IWerner etall 
I2OO7I : ISanders fc Fab ian"2007'V Cosmic rays have also been 
inferred from the excess abundance of ®Li in metal-poor halo 
stars, since ^ Li could be produ ced in spallation reactions by 
cosmic rays (|Nath et al.ll2006l '). 

The cosmic-ray heating of t he ICM has been 
studied by several authors (iBoehringer fc Morfilll 



19881: iLoewenstein et all Il99ll : iRephaeU fc Silkl 



Colafrancesco et al.l 2004 : lJubelgas et al 



1995 



2006 



Pfrommer et al.ll2006l '). While many studies have found that 



cosmic rays could in principle be dynamically important, 
none have constructed successful models in which CR 
heating prevents a cooling catastroph e. In the steady-state 
model of IBoehringer fc Morfilll (Il988l), a significant c ooling 
fiow (~ SOOMQyr"^) developed. Loewenstein et al.l (|l99ll ) 
proposed a hydrostatic model where radiative cooling is 
fully balanced by hydromagnetic-wave-mediated cosmic-ray 
heating and thermal conduction. Their model could not 
fit the observational gas density profile; in particular, 
they found that the intra-cluster medium would quickly 
become CR pressure dominated at a level inconsistent with 
observations, lon g befor e heating effects become significant. 
iRephaeU fc sil^ (|l995h estimated the CR heating rate 
from Coulomb collisions alone, which they argued could 
be significant, but did not construct a specific model for 
the intraclust er gas, which could be co mpared against 
observations. IColafrancesco et al.l (|2004h do construct 
semi-analytic models of cluster density and temperature 
profiles, which differ significantly from ours: they only 
consider Coulomb and hadronic heating; the ICM is not 
in thermal equilibrium but evolves strongly as a function 
of timsQ; 

they do not solve for hydrostatic equilibrium, 
and thus their input density profiles ne(r) do not evolve 
as the temperature profi le T(r) evolves. More recently, 
lJubelgas et al.l (|2006l) and iPfrommer et all (|2006l) have run 
a suite of high-resolution 3D numerical simulations analyz- 
ing the role of cosmic rays in clusters; these simulations 
represent the most comprehensive and careful treatment 
of this problem to date. These authors find that cosmic 
ray heating cannot stem a cooling fiow; in particular they 
find that the increased compressibility (due to the softer 
adiabatic index of CRs) can lead to enhanced cooling. 
However, as they note, this may be due to the relatively 
modest production of CRs in their self-consistent treatment, 
where cooling gas gives rise to star formation and hence 
supernovae (the source of CRs in their model). Moreover, 
they only consider cosmic ray heating through Coulomb 
interactions with the ICM, which is much less than the 
hydromagnetic-wave mediated cosmic-ray heating in our 
models. The inclusion of AGN-mediated CR production 



^ This may be hard to reconcile with observational evidence that 
the temperature profiles of the intracluster gas are well described 
by an approxima tely universal ma thematical function across a 
range in redshift l l Allen et al.ll200ll') . 
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and hydromagnetic wave mediated cosmic-ray heating 
could significantly alter this picture. 

In contrast to previous work, we find that it is indeed 
possible to construct models in which cosmic-ray heating 
strongly suppresses the cooling flow. Although our results 
are generic, we choose the typical cool core cluster Abell 
2199 as our fiducial model. By using ID numerical hydro- 
dynamic simulations, we demonstrate that, starting from 
a state far from thermal equilibrium (isothermal tempera- 
ture profile), the cluster will relax to a stable steady quasi- 
equilibrium state, in which the accretion rate is highly re- 
duced and the temperature and density profiles are con- 
sistent with observations. We take the relevant cosmic ray 
physics into account: various CR energy loss mechanisms, 
heating of the ICM by cosmic rays, cosmic ray pressure sup- 
port, and cosmic ray transport (advection and diffusion, see 
§ 12. 2|) . We also incorporate a moderate level of thermal con- 
duction into our models. As our simulations show, thermal 
conduction delays the onset of the cooling catastrophe dur- 
ing the early stages of the ICM evolution, while the feedback 
heating by cosmic rays suppresses the cooling catastrophe 
quickly after it starts. 

The remainder of this paper is organized as follows. In 
§ [2] we describe cosmic-ray heating mechanisms and some 
physical assumptions we made in our models. In §|3]we write 
down the relevant differential equations of the intracluster 
medium. In § 2] we present the results of a series of numeri- 
cal hydrodynamic simulations that investigate the evolution 
of the ICM, starting from isothermal hydrostatic equilib- 
rium. We conclude with a summary of our results in § [S] 
The derivation of the time-dependent cosmic-ray equations 
is presented in the Appendix. In computing luminosity and 
angular diameter distances, we have rescaled observational 
results if the original paper used a different cosmology to 
the values f2m = 0.3, Ha = 0.7, h = 0.7. 



2 COSMIC RAY PHYSICS IN CLUSTERS OF 
GALAXIES 

2.1 Cosmic-ray heating mechanisms 

In this subsection, we quantify various cosmic ray energy loss 
mechanisms in clusters of galaxies, and their heating effects. 
For simplicity, we only consider relativistic protons in cosmic 
rays. Relativistic electrons are generally insufficient to heat 
the intracluster gas (|Rephaell.l987i ). 



2.1.1 Hydromagnetic waves 

The significance of the resonant interaction of cosmic rays 
with hydromagnetic waves has been recognized and dis- 
cusse d by ma n y au t hors in various astrophysical con texts 
fe.g.,lMel rosc fl9 68|); ISkilliii^ (Il97ll): Blandford & Ostriker' 
1 19781 '): ILoewenstein et all (|l99ll ): iBrunetti et al.. (,2004) 1. 
During wave-particle resonance, the waves may be damped 
or grow exponentially, dependin g on t he cosmic r ay dis - 
tribution function (see iMehosd l|l968l '): iKulsrudI (|2005h . 
Chap. 12). If the distribution function of cosmic rays 
is sufficiently anisotropic due to their streaming motion 

i drive n by the cosmic-ray pressure gradient; see iKulsrudI 
20051 ) ■ Chap. 12.6) with a streaming speed in excess of 



the local Alfven speed, they will render "forward" Alfven 
waves propagating along the magnetic lines in the direc- 
tion of the strea ming unstable, while backward Alfve n waves 
are d a mped ( | Kulsrudl (i2005i'). C hap. 12; see also iLerchd 
(|l967^ ■ [Kuisrud fc Feared (Il969l )'). During this quasi-linear 
cyclotron resonance (here also called the cosmic ray stream- 
ing instability), f orward (and nearly forwa rd) Alfven waves 
grow fastest (e.g. iKulsrud fc Feared ri969l ). scatter the cos- 
mic rays and reduce the cos mic-ray stream ing speed to 
around the Alfven speed (e.g. ISkillingj|l971^ . Cosmic rays 
are thus advected with these waves as scattering centers. 

These unstable Alfven waves grow exponentially un- 
til saturated by nonlinear processes , e.g., nonlinear Lan- 



dau d a mping (see Lee fc Volkl (Il973|): ICesarskv fc KulsrudI 
(|l98ll ): rVoelk et al.l (| 1984 ): IKulsrudI (|2005l ). Chap. 11.5) and 

are thus dissipated efficiently in the ionized thermal gas. 
Note that the decay of a forward Alfven wave to a back- 
ward Alfven wave and a forward sound wave is forbidden 
in the ICM where the sound speed Ca is usually larger than 
the Alfven spee d va (the requi rement for the wa ve decay 
is DA > Cs; see IKulsrudI (|2005l ). Chap. 11.4 and ISkiUind 
(19751^)). Therefore, momentum and energy are transferred 
from relativistic protons to the waves and hence to the 
intracluster medium at the rates (per volume) |VPc| and 
— (u + Va) -V Pc, resp e ctivel y (see Appendix A for the de- 
tails; also see IWentzell (|l97ll )). where Pc is the cosmic-ray 
pressure, u is the bulk velocity of the thermal plasma and 
Da is the propagation velocity of the hydromagnetic waves 
relative to the plasma, which is equal to the local Alfven 
velocity of the gas. The energy gai ned from the co smic rays 
can accelerate and heat the gas (IWentze]| [l97ll ). The in- 
crease in gas kinetic energy due to the work done by the CR 
pressure gradient is obviously —u -V Pc. Thus the cosmic- 
ray heating rate of the ICM due to the dissipation of hy- 
dromagnetic wav es (hereinafter designated as "cosmic-ra : 
wave heating") is (| Wentzel|[l97ll : iBoehringer fc Morfiii|[l98; 



ILoewenstein et~al 



19911 ) 



= -VA 



■VPc 



(1) 



We note that only Alfven waves self-excited by the 
cosmic-ray streaming instability are considered in this pa- 
per. In reality, the MHD waves and turbulence in galaxy 
clusters may be much more complex. The MHD turbulence 
driven by cluster mergers may reaccelerate cosmic rays, 
which has been studied by many authors to explain radio 
halos and hard X-ray tails in some galaxy clusters (e.g. 
FetrosianI 1200 ll : IBrunetti et~al] 12004 IBrunetti fc LazarianI 



20071 ). These models usually assume an isotropic, homoge- 



neous cosmic-ray phase space distribution function fp{p,t), 
and the cosmic rays are generally reaccelerated through the 
wave-particle resonance as long as dfji/d p < (the wav es 
grow exponentially if dfp/dp > 0; see, e.g.. lMelrosd (|l968l )). 
In our model, the cosmic rays are mainly injected into the 
ICM by the central AGN and the cosmic ray pressure gradi- 
ent is strong. In this case, the cosmic-ray streaming along the 
CR pressure gradient is significant, and the forward Alfven 
waves self-excited by the cosmic-ray streaming instability 
should be the main MHD waves res ponsible for the cosmic- 
ray scattering (see IChandranI |2000| ) . Note that cosmic ray 
streaming also depends on the details of the CR scattering 
off the small-scale MHD turbulence in the ICM. The details 
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of the latter is still poorly understood at this point, and 
further observational and theoretical work is needed. 

It is worth noting, not only the "drifted" cosmic- 
ray distributions, but also prolate (pz > p±) or oblate 
(pz < p±) CR distributions may generate Alfve n waves 
throu gh cyclotron resonance (see iKulsrudI (|200d ). Chap. 
12.3). lLazarian fc BeresnvakI l|2006l ) recently proposed a new 
model for the generation of small-scale Alfven waves through 
cyclotron resonance, where the prolate or oblate distribu- 
tions of cosmic rays are driven by turbulent compressions of 
magnetic field. 



2.1.2 Coulomb interactions 

The relativistic protons can transfer energy to the gas via 
Coulomb collisions with the ambient ionized gas. The heat- 
ing rate due to Coulomb interactions of a fully ionized gas 
with a cosmic-ray particle with c harge Ze, velocity v = l3c 
and kinetic energy E is give by (jMannheim fc Schlickeisen 
1 19941 ) 



^] =4.96x10- 



\cm~3 / 



13' 



cm~3 / + X- 



-ergs s 



(2) 



where = 0.0286(r/2 x 10''K)^''^ T and are the ambi- 
ent electron temperature and number density. For simplicity, 
we neglect Xm in our calculations, which is a good approxi- 
mation for /3 > 0.15 (i.e. E > 10 MeV) for a typical cluster 
temperature T ~ 5 keV. 



2.1.3 Hadronic collisions 

The cosmic-ray protons interact hadronically with the am- 
bient thermal gas and produce mainly tt"*", tv~ , vr", pro- 
vided their kinetic energy exceeds the threshold Ethi ~ 282 
MeV for the reaction. The neutral pious decay after a 
mean lifetime of 9 x lO"^'^ s into 7-rays (tt" -> 27), which 
may be detected by 7-ray observations with imaging at- 
mospheric Cherenkov telescopes or the GLAST satellitfl 
One may thus test our model with the future observation of 
these 7r°-decay induced 7- rays from galaxy clusters (e.g., see 
iHinton fc Domainkd[2007h . The charged pions decay into e 
and neutrinos (7r='= fi"^ + ly^/V^ e"^ + v^jva + f m + '^i^)- 
Since the limiting value of the inelasticity of these hadroni c 
collisions is roughly 1/2 (jMannheim fc Schlickeiseil Il994h . 
the energy loss rate of a c osmic-ray proton d ue to pion pro- 
duction is approximately l|Enfilin et al.|[2006l ) 



~dt 



-0.5nKapp/3cEe{E ~ Etbr 



(3) 



where riN = no/(l — ^Y) is the target nucleon density in 
the ICM, Y is the helium fraction, and o-pp is the pp cross 
section for the incident proton. We adopt an approximate 
value for o-pp (gpp ^ 37.2 mbarn) by using equation (69) of 
lEnfihn et all (|2006l ) and assuming that the spectral index of 
the TT^-decay induced 7-ray spectrum is 2.5. 

To estimate the total energy loss rate of CRs due to 



^ Gamma-ray Large Array Space Telescope (GLAST), homepage 
|http : / /glast . gsfc . nasa. govT] 



Coulomb interactions and hadronic collisions, we need to de- 
termine the cosmic-ray energy spectrum. Galactic CR obser- 
vations and many CR acceleration models suggest that the 
CR s pectrum is a power-law in momentum (see Schlickeiseil 
|2002| . for a review). However, in low energy regimes, the en- 
ergy losses of cosmic rays are dominated by Coulomb inter- 
actions, which subs tantially flatten the spectrum. Similar to 
lEnfihn et al.1 l|2006l ). we derive an approximate steady-state 
CR spectrum to estimate the total Coulomb and hadronic 
loss rates. Assuming that the cosmic rays are injected con- 
tinuously and subject to energy losses through Coulomb and 
hadronic collisions, the cosmic ray spectrum in a homoge- 
neous environment follows the evolution equation: 



dfE{E,t) d ( dE 



Q^{E) , (4) 



where the cosmic-ray spectrum fE{E,t) is defined as 
fE{E,t) = 

and 



dN . 2 , / X dp 
^=4V/.(.,P,t)-^, (5) 



dE 
It 



dE 
Ht 



dE 
'dt 



(6) 

c \ / h 

is the energy loss rate due to Coulomb and hadronic colli- 
sions. For simplicity, here the cosmic-ray phase space dis- 
tribution function fp{x,p,t) is assumed to be isotropic in 
momentum space and we approximate hadronic losses as 
continuous, rather than impulsive. Assuming that the CR 
injection spectrum is a power law in momentum, i.e. 



]p{p) ocp °'e{p-pi), 



(7) 



where ^(a;) is the Heaviside step function and pi is the 
momentum cutoff, we can get the CR injection spectrum 
QTi{E): 

Qe{E) <x{E + Eo)[E{E + 2Eo)]-^"+'^/'e{E - Ei) , (8) 

where Ei is the energy cutofi^ and Eo — 938 MeV is the pro- 
ton rest energy. We find the asymptotic steady-state spec- 
trum by assuming negligible hadronic and Coulomb losses 
in the low and high energy regimes, respectively: 



for E > E,, 

for El <E <^ E,, 



(9) 



^""^^^-^"[{E/E^r^-^y^ 

where Acr is the normalization factor which can be deter- 
mined in terms of the energy density (i?c) of cosmic rays, 
and the cross-over energy E^ « 706 MeV depends on the 
ratio of the Coulomb to hadronic loss rates. We note that 
a similar asymptotic statio nary C R sp ectrum has al s o bee n 
derived by iBrunetti et all (|2004h and lEnfilin et aP l|2006h . 
We adopt a simple analytic approximation for the steady- 
state CR spectrum with the same asymptotic behaviors: 

f (F\ = AcrOjE - El) 

■""^ ' (£/£;.)" (£;/£.)("-2)/2 ■ ^ ' 

Using the approximate steady-state spectrum (eq. |10)). 
we get the overall Coulomb loss rate of cosmic rays 



Tc 



dE 



-1.65 x 10" 



(^) 

V cm / 



Ec 



ergs cm" 



ergs s cm 



(11) 
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where Ec is the energy density in cosmic rays, and we adopt 
a = 2.5 and Ei = 10 MeV. The loss rate Tc depends slightly 
on the value of Ei , but not very sensitively, since in the low 
energy regimes fsi^E) is quite flat (i.e., the CR spectral in- 
dex, (q — 2)/2, is quite small). Similarly, the overall hadronic 
loss rate per volume is 



-5.86 X 10" 



16 / 



ergs cm 



ergs s ^ cm ^ 



(12) 



We note that CR Coulomb losses are usually sub-dominant 
with respect to hadronic losses if we use the stationary CR 
spectrum (eq. [10]) with a < 3.7 to calculate the CR energy 
loss rates. For higher values of a or using the CR injection 
spectrum (eq. [7] with a > 2) instead of the stationary spec- 
trum. Coulomb losses usually dominate. 

Therefore, the total cosmic-ray energy loss rate due 
to Coulomb and hadronic collisions is Fioss = Th -I- Fc = 
—(cneEc where ~ 7.51 x 10~^^ cm^ s~^ is the rate coef- 
ficient for coUisional energy loss of the cosmic rays. While 
the CR energy lost in Coulomb interactions heats the ICM, 
most of the CR energy lost in hadronic collisions tends to 
escape via gamma rays and neutri nos. During hadronic colli- 
sions, a small fraction {^^ 1/6, see lMannheim fc Schlickeiserl 
of the inelastic energy goes into secondary electrons. 
While the cooling of high-energy electrons (7 > 10'') is dom- 
inated by synchrotron and inverse Compton losses, most 
of the mildly relativistic (7 < 200) electrons will heat the 
ICM by Coulomb interactions with the plasma electrons and 
thro ugh plasma osc illations and excitation of Alfven waves 
(e.g. lRephaeli|[l979l ). Since the spectrum of secondary elec- 
trons is dominated by mildly relativistic electrons, here we 
assume that these secondary electrons lose most of their en- 
ergy through thermalization and thus heat the ICM. There- 
fore, the ICM is heated by cosmic rays through Coulomb 
and hadronic collisions (hereinafter designated as "CR col- 
lisional heating") at a rate of Fcoii ~ — Fc — Fii/6 = r/cJicSc, 
where 77^ = 2.63 x 10"^*^ cm^ s"\ 

Note that the real CR spectrum in galaxy clusters may 
differ substantially with our simple steady state spectrum, 
since our calculation neglects the influence of CR transport 
and energy losses due to the generation of hydromagnetic 
waves. The detailed shape of the spectrum has relatively 
little impact on hadronic losses. The Coulomb heating rate 
does depend on the CR spectrum, but even if these effects 
are included, the CR spectrum is still strongly flattened 
at low energies due to strong Coulomb losses. The heating 
rate from the steady-state CR spectrum (which is quickly 
reached since tcouiomb <C ^duster at low energies) is an order 
of magnitude less than the heating rate calculated from t he 
injection spectrum (e.g., equation (3) in lNath et al.l l|2006l) ). 
which is only a short transient. As shown in § 14.41 CR col- 



lisional heating (?7c = 2.63 x 10~ cm s~ ), is negligible 
comparing to CR wave heating. Thus, our results should 
not change appreciably if a new coUisional heating rate cal- 
culated with a more realistic CR spectrum is used instead. 



2.2 Cosmic-ray propagation and assumptions 

In our models, cosmic rays are assumed to be primarily in- 
jected into the ICM by AGN-produced radio cocoons or bub- 
bles. Cosmic rays then propagate along the magnetic field 
lines, which, for simplicity, are assumed to be mainly radial 
on a large scale. Such a magnetic field structure could be 
created by the bubbles themselves: magnetic fields stretch 
out behind the bubble and becom e approximately paralle l 
to the (radial) direction of motion (iRuszk owski et ai1l2007l ). 
Cosmic rays can thus be viewed as streaming radially on a 
large scale. If so, cosmic ray propagation can be described 
approximately by a simplified one-dimensional model with 
spherical symmetry. In the future, it will be interesting to 
consider more realistic 3D hydromagnetic simulations, which 
(among other effects) could potentially take into account the 
effects of cross-field diffusion of CRs. 

The time-dependent cosmic-ray equations governing 
CR transport and CR energy loss due to the generation of 
Alfven waves are derived in Appendix A. As shown in the 
previous subsection, cosmic rays in the ICM can also lose 
their energy through Coulomb and hadronic collisions with 
the ambient thermal gas, at a rate Fioss = (cneEc- Taking 
CR coUisional losses into account, the net CR source func- 
tion Q in equation HA15|I may be written as Q = Qc— Cc^ieSc, 
where Qc is the source (injection) function of cosmic-ray en- 
ergy. 

Assuming spherical symmetry, the CR energy equation 
(|A15|) may be rewritten as 



dE, 
dt 



(7c-1)(u+i'a) 



dEc 
dr 



-CcTlcEc 



1 djr^F,) 

J.2 



+Qc, (13) 



where u is the radial velocity of the thermal plasma, «a is 
the local Alfven speed of the thermal gas, is the cosmic- 
ray energy Hux, and 7c = 4/3 is the adiabatic index for the 
cosmic rays. The cosmic-ray energy flux Fc is (see eq. [A14j l 



Fc = 7c(?i + va)Ec ■ 



dEc 
dr ' 



(14) 



where Kc is the diffusion coefficient of cosmic rays. The two 
terms in equation (|14p clearly represent the spread of cosmic 
rays by advection and diffusion respectively. 

Magnetic fields in the ICM have been measured us- 
ing a variety of techniques, such as Faraday rotation mea- 
surements and studies of sync hrotron radiation and inverse 
Compton X-ray emission (see ICarilli fc Tavloil [2002 . for a 
review). These measurements imply that the ICM of most 
clusters is substantially magnetized, with a typical field 
strength of order 1 /xG with high areal filling factors out to 
Mpc radii. In the cores of cool core clusters, these measure- 
ments also suggest that magnetic f ield strength is ty pically 
much higher, up to 10s of ^G fe.g. [Alien et al.ll2001 ). 

The magnetic field is usually not uniform and its distri- 
bution in galax y clusters is still far fr om clear. For definite- 
ness, similar to [jubelgas et al.l (|2006l ). we assume that the 
magnetic pressure (Pb = B^/{8tt)) is a fixed fraction of the 
thermal pressure (Pg), which corresponds to 



(15) 



In the rest of this paper, we assume that Ps/Pg ~ 0.06, 
which corresponds to B ~ 10 — 20 /xG in the central re- 
gions of the cluster in the final steady state. Such magnetic 
field amplitudes have been found in the cores of cool core 
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clusters through the analysis of Farad ay Rotation measure 
(RM) maps (e.g. IVogt fc Enfiiiiill2003h . We neglect pressure 
support from magnetic fields, since this small level of mag- 
netic field pressure will not substantially change our main 
results. We also neglect the pressure of the hydromagnetic 
waves, which should at most be comparable to the pressure 
of the main field. 

The diffusion coefficient of cosmic rays in the intra- 
cluster medium, which depends on the frequecy of pitch- 
angle scattering and the cosmic-ray momentum spectrum, 
is highly uncertain and may vary substantially between dif- 
ferent parts of galaxy clusters. F or definiteness, here w e fol- 
low a simple treatment given bv lJubelgas et al] l|2006l ) and 
discuss the dependence of our model on radial profiles of Kc 
in § 14.41 In Kolmogorov-like MFI D turbulence , the parallel 
diffusivity is expected to scale as l|Enfilinll2003l ) 



2/3 D-1/3 



(16) 



where Ib gives a characteristic length scale for the magnetic 
field of strength B. If ere we have ignored the weak depen- 
dence of Kc on the cosmic-ray energy distribution. Ib is fairly 
unknown in galaxy clusters. For defini teness, we assume tha t 
Ib scales with the local Jeans length jjubelgas et al.|[200d ), 
and thus obtain: 



0.12 cm- 



-1/2 



T 



1.6 keV 



1/6 



'\ (17) 



where we fixed the normalization kq by assuming that 
Kc at the center of the fiducial cluster Abell 2199 equals 
to the estimated diffusion coefficient along the magnetic 
field lines in the interstellar medi um of our own Galaxy : 
«o = Kc, ISM « 3 X 10^* cm^ s"^ l|Berezinskii et aLllTooa ). 
The values of ne(rinin) = 0.12 cm~^ and T(rmin) = 1.6 keV 
at the inner boundary rmin ~ 1 kpc are obtained di- 
rectly from Chandra o bservations of the cluster Abell 2199 
ij Johnstone et al.]l2002l l. Obviously, this estimate is highly 
uncertain, but it is consistent with the approximate bound 
of ^ 7.5 X lO^Viib^^st-^^ 7 cm s required for large 
cavities of radii riobe ~ 5nobc,5 kpc to remain stable to the 
diffusion of cosmic ray s in time tiobe ^ lO^tiobe.7 y r, as re- 
quired by observations jMathews fc Brighenti|[2007^ , as well 
as theoretical estimates l|Enfilinll2003 v 



2.3 Injection of cosmic rays by a central AGN 

The injection of CRs is specified by the cosmic-ray source 
function Qc, which depends on the spatial profile of CR in- 
jection from AGN activity, supernovae, structure formation 
shocks and merger shocks. The form of Qc is far from clear. 
In this paper, we only consider CR injection from the cen- 
tral AGN in clusters, via jet-ICM interactions, which create 
buoyant bubbles filled with cosmic rays. When the bubbles 
rise, they inject cosmic rays into the ICM through diffu- 
sion or the shredding of the bubbles by Rayleigh- Taylor 
(RT) and Kelvin-Helmholtz (KH) instabilities. This is an 
efficient means of transporting cosmic rays from the AGN 
out to large distances: the buoyancy timescale is typically 
comparable to (at most several times) the sound crossing 
time tsc ~ 10*riooc~loooyi' ^ radius r ~ lOOrioo kpc 
and sound speed Ca ~ 1000cs,iooo kms~^ (e-g-, see table 3 in 
iBirzan et al.l l|2004l )). whereas the CR diffusion time along 



field lines is idifiusion ~ 3 x lO^^^r^QQAt;. 283"" (where the diffu- 
sivity Kc ^ 10^®Kc,28 cm^ s^^), too loug to affect the ther- 
mal state of gas at the cooling radius. Efficient advection of 
cosmic rays decreases this timescale by about an order of 
magnitude (see tj4.4|l . still too slow. Bubble transport is a 
key ingredient of our model: if excluded, cosmic-ray trans- 
port is too slow to allow sig nificant heating (e.g., models of 
iBoehringer fc Morfilll (|l98^ )'). 

AGN activity is likely to be intermittent on a timescale 
of order the Salpeter time ts ^ 10^ yr, and pos s ibly a s short 
as ~ 10** - lO'^ yr (Ruszkowski & Begelman| (|2002l '). here- 
after RB02; [Reynold s & Bcgelman (1997)), which is much 
shorter than the bubble rising time. Note that the bubble 
rise time is usually much shorter than the gas cooling time 
(RB02). It is thus justifiable to assume that the CR injection 
into the ICM from the buoyant bubbles, which are produced 
by a succession of AGN outbursts, can be treated in a time- 
averaged sense. Since CR transport timescales are much 
shorter than thermal timescales, we assume that the cosmic 
rays are injected into the ICM instantaneously and neglect 
any delay between central AGN activity and the cosmic-ray 
injection (similar to instantaneous AGN mechani cal heating 
models, e.g., RB02. iBrighenti fc Mathews! (120031 )'). 

The rate at which bubbles are disrupted is highly 
uncertain, since the nature of the physical mechanism 
which protects them from RT and KH instabilities — 
perhaps an ordered magnetic field at the bubble sur- 
face, or thermal conduction /plasma viscosity (|Kaiser et al.l 
l2005l : [Reynolds etall l2005h . or even the initial decelera- 
tion and drag on the bub ble-ICM interface during inflation 
(|Pizzolato fc Sokeilliooel ') -is not well understood. Likewise, 
the diffusion rate of CRs out of the bubbles is highly uncer- 
tain, particularly given the unk nown magnetic field to pol- 
ogy at the bubble interface (see iRuszkowski et al.ll2007l , for 
recent simulations). Nonetheless, the observational fact re- 
mains that bubbles are seen to survive intact to large radii 
(with an av erage projected rad ius of ~ 20 kpc in a sample of 
16 clusters (|Birzan et al.|[20o3 )). We hence parametrize our 
ignorance of the bubble disruption rate by simply assum- 
ing that the spherically integrated CR energy flux in the 
buoyant bubbles is a power law with radius: 



-'^buh 



for r > ro. 



(18) 



where e is the efficiency with which the rest-mass energy 
of the ICM cooling flow is converted into the cosmic-ray 
energy in the bubbles, Min is the mass accretion rate across 
the inner boundary of the simulation, c is the speed of light, 
!^ is a positive constant, and ro is a characteristic radius 
where the bubbles are created. The decline of I/bubbio with 
radius reflects the transfer of CRs to the ICM. In particular, 
the cosmic-ray energy injection rate into the ICM per unit 
volume is given by: 

„ „ 1 9I/bubblc 

Qc = -V-Fbubbio ~ ~- 1 1 - e 



4j^j.2 

ro 



47rr|j 



dr 

-3- 



1-e- 



-{r/ro)' 



(19) 



where Fbubbic ~ I'bubbiG/(47rr^), and we have introduced an 
inner injection cutoff term, which reflects the finite radius 
r ~ ro at which bubbles are injected; it is very similar to the 
AGN heating cutoff term in RB02. Here ro was taken to be 
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20 kpc. Inclusion of this term allows us to apply equation 
(|19p at all radii in our simulations. 

Besides CR injection into the ICM, bubbles also lose 
energy by expansion as they rise. Therefore, strictly speak- 
ing, Lbubbio defined in equation (I18|) is only a means of 
parametrizing the spatial distribution of the CR injection 
rate (eq. dH), and is generally different from the total en- 
ergy flux of cosmic rays. We can include the PdV work from 
bubble expansion in equation ((THJ: 



Lbui 



pdv 



for r > ro. (20) 



As before, the first term on the right-hand side of equation 
H20p represents the cosmic-ray injection into the ICM, while 
the second term jFpdv represents the PdV work done on the 
ICM during the bubble expansion. The AGN feedback effi- 
ciency etot is: 



etot 



□ (ro) 



A/i„c2 



(21) 



If the bubble is nearly in pressure equilibrium with the ICM, 
.T^pdv satisfies the equation: 



(7. 



e(r/ro)"" + J^pdv 



pdv 



dr 



(22) 



If we further assume that Pg scales approximately with r 
(r > 0), the solution of equation (|22|l is: 



pdv ' 



(ttot 



-T I UJ 



+ a(r) 



where uj = "jc/isic — 1), and 



LOV — T 



tT(^-l) 



— t/uj 



In 



if iy/ i 



if 1/ 



(23) 



(24) 



The AGN mechanical heating model of RB02 assumes 
no cosmic ray leakage from the bubbles and finds that J-pdv 
scales with Pg''" ^^^''^ ^ which is the solution of equation (|22p 
when e = 0. On the other hand, in this paper, we consider 
the heating of the ICM by the cosmic rays leaked from these 
AGN-produced bubbles, and assume that e ^ etot- In this 
case, the effect of the bubble expansion or contraction is de- 
scribed by the variation of Q{r) with radius. For < < t, 
Q(r) first decreases and then increases with radius, which 
means that the bubble expands first and contracts later. For 
V > T, the bubble contracts first and expands later. We are 
not interested in the models with high values oi v {v >^ 1.0), 
where the cosmic rays are essentially dispersed into the ICM 
at the cluster center. When v is very small (u < 0.1), the 
CR injection is less important and the bubble PdV work 
(the model of RB02) dominates. In this paper, we are in- 
terested in the models with a moderate level of the bubble 
disruption (0.1 ^ v <^ 1.0), where the CR injection into the 
ICM dominates and hence the CR heating of the ICM may 
be significant. Figure [T] shows the ratio of the bubble PdV 
work rate to the CR injection rate for different values of 
y and r (t ~ 0.5 — 1.0 in the central regions of the clus- 
ter Abell 2199). As can be clearly seen, bubble expansion is 
subdominant for u — 0.3 or 0.7, and becomes comparable to 
the cosmic ray injection only for very fiat disruption profiles. 




Figure 1. The ratio of the bubble PdV work rate to the CR 
injection rate, plotted as a function of radius (7c is taken to be 
4/3). For each line type, the upper line (green) corresponds to 
T = 1, while the lower one (black) corresponds to r = 0.5. Note 
that, in the central regions of the cluster, the bubble expansion 
is less important for smaller values of t. A positive value of the 
PdV work rate means that the bubble expands, while a negative 
value corresponds to the bubble contraction. 



v = 0.1, reverting to the case studied by RB02. Hence, we 
shall neglect bubble expansion in this paper. 

The cosmic ray injection described by equations (|18|l 
and H19|l is obviously simplified, and would be worth refin- 
ing in more detail once more is known both observationally 
and theoretically about bubble disruption. Nonetheless, our 
model should be relatively robust to the details of the bub- 
ble disruption profile. For instance, as we show in 14.41 our 
model is very robust to the value of the model parameters 
(i.e., no fine-tuning of the model parameters e, v needed). 



3 BASIC EQUATIONS 

In our spherically symmetric model of galaxy clusters, the 
dark matter dis tribution is given by a Navarro - Frenk- White 
(NFW) profile jNavarro. Frenk fc Whitj|l997l ): 



PDM(r-) 



Mo/2tv 
rir -\- rsY' 



(25) 



where Vs is the standard scale radius of the NFW profile and 
Afo is a characteristic mass. 

In this paper, we will take the cooling flow cluster 
Abell 2199 as our flducial model to study cosmic-ray heat- 
ing in galaxy clusters. For this cluster, the parameters of 
the NFW proflle are Mp = 3.8 x 10^"* Mq, = 390 kpc 
(|Zakamska fc Naravan|[2003l l. Since we are interested in the 
cluster within the cooling radius, the central cD galaxy 
NGC 6166 may be dynamically important as well. We adopt 
a King profile with core radius and one-dimensional ve- 
locity dispersion of rg = 2.83 kpc and a = 200 km s~^, 
resp ectively, for t he density distribution of NGC 6166 
(Kelson et al.ll2002l ): 

PcD(r) = 



[14- (r/rg)2]3/2' 
where po is the central density. 



9a' 



(26) 



(27) 
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The intracluster medium is subject to radiative cooling, 
thermal conduction and cosmic-ray heating. The governing 
hydrodynamic equations are 

dp 
dt 



+ V-(pu) 







(28) 



— + V-(5u) = -VPg- VPc-pV$ , (29) 
ot 



at 



+ V-(£gu) 



Pg V ■ u - V • F - ntA{T) + rjcUcEc 



(30) 



where p is the gas density, Pg is the gas pressure, Eg is the 
gas energy density, S — pu is the gas momentum vector, 
Pc = (7c — l)Ec is the cosmic-ray pressure and "I> is the 
gravitational potential, which is obtained by solving Pois- 
son's equation: 



47rG(pE 



' PcDj 



(31) 



where we have neglected the self-gravity of the intracluster 
medium (ICM). Hence, ^ can be written as $ = "1?dm -I-'J'cD, 
where 



$1 



2GAfoln(l -hr/re 



rs r/rs 
is the cluster dark matter potential and 



ln(r/rg + + (r/rg)2) 
r/rg 



(32) 



(33) 



is the gravitational potential contributed by the central cD 
galaxy. In our spherically symmetric model, we define the 
gravitation acceleration g = d^/dr. 
We adopt the ideal gas law, 

pksT 



P. = 



= ncKBl , 



(34) 

pm^ p 

where fca is Boltzmann's constant, is the atomic mass 
unit, and p and pc are the mean molecular weight per ther- 
mal particle and per electron, respectively. We assume that 
the gas is fully ionized with hydrogen fracti on X = 0. 7 
and helium fraction Y = 0.28 (fZakamska fc Naravanll2003l ). 
so that p, = 0.62 and Pe ~ 1.18. We use an analytic fit 
iTozzi fc Norman 2001) to the cooling fun ction based on 
calculations by Sutherland fc Dopital l|l993l ). 



n^A(r) =1.0 X 10"^^ 



\ cm ^ / \ cm / 



-3 -1 

ergs cm s 



(35) 



where is the ion number density. For an average metal- 
licity Z — 0.3^0, the constants are 5i = —1.7, S2 = 0.5, 
Ci = 8.6 X 10-^ C2 = 5.8 X 10"^ and C3 = 6.3 x 10-^ 
and we can approximate niric — 0.704(p/mp)^, where rrip is 
the proton mass. We manually truncate the cooling below 
a minimum temperature of 0.0 3 keV, since equation II35P is 
only vahd for ksT > 0.03 keV l|Tozzi fc NormanlbOOll ). 

In equation (|30[) . F is the heat flux due to electron 
conduction, 

F = -fuspVT, (36) 
where / (0 ^ / ^ 1) is a conductivity reduction factor due 



Table 1. List of Simulations. 



Run Heating 



fl 



A 


None 




N/A 


N/A 


N/A 


Bl 


Conduction 




0.4 


N/A 


N/A 


B2 


Conduction 




0.8 


N/A 


N/A 


B3 


Conduction 




0.6 


N/A 


N/A 


C 


Conduction, 


CRf, 


0.3 


0.003 


0.3 


Dl 


Conduction, 


CR~ 


0.1 


0.003 


0.3 


D2 


Conduction, 


CR 


0.4 


0.003 


0.3 


El 


Conduction, 


CR 


0.3 


0.05 


0.3 


E2 


Conduction, 


CR 


0.3 


0.0003 


0.3 


Fl 


Conduction, 


CR 


0.3 


0.003 


0.1 


F2 


Conduction, 


CR 


0.3 


0.003 


0.7 


F3 


Conduction, 


CR 


0.3 


0.003 


1.5 


Gl^ 

G2 / 


Conduction, 


CR 


0.3 


0.003 


0.3 


Conduction, 


CR 


0.3 


0.003 


0.3 


G3£ 


Conduction, 


CR 


0.3 


0.003 


0.3 













Conductivity suppression factor relative to the Spitzer 
value. 

Efficiency of cosmic ray injection due to accretion-triggered 
AGN activity. See eqs. JTS)! and ITOt . 
^ See eqs. JS} and ((T9j. 

Cosmic-ray heating. 
^ Runs Gl, G2 and G3 are performed to check the dependence 
of our model on the CR diffusion coefficient. For run Gl, Kc 
has the form of equation 1171 1 with kq = 3 x lO'^^ cm^ s~^. 
f For run G2, Kc has the form of equation I I17II with kq = 
3 X 10^^ cm^ s-l. 

9 For run G3, Kc is taken to be constant throughout the clus- 
ter: Kc = 3 X 10^® cm^ s~i. 



to magnetic fi eld suppressio n and nsp is the classical Spitzer 



to magnetic n eia suppressio r 
conductivity (|Spitzerlll962h . 



1.84 X 10" 
iiiA 



rn5/2 -1 -1 

-1 ' ergs s K ' cm , 



(37) 



where InA ~ 37 is the usual Coulomb logarithm. For sim- 
plicity, here we assume that / is constant throughout the 
cluster. In real clusters, heat transport may be much more 
complex, depending on the plasma magnetization and tur- 
bulence driving. Heat transport through turbulent fluid mo- 
tions may also need to be taken into account (see iLazarianl 
I2OO6I ). Since at present there is no consensus on the nature 
of conductivity in a turbulent magnetized plasma, we adopt 
the same assumption of Spitzer conductivity (with a con- 
stant suppression factor) that most authors do. We do show 
that (unlike others) no flne-tuning of / is necessary in our 
models. 



4 HYDRODYNAMIC SIMULATIONS 

In this section, we follow the long-term evolution of a series 
of spherically symmetric cluster models from a state far from 
equilibrium to investigate whether the cluster will relax to 
a stable quasi-equilibrium state. Basic information for the 
set of simulations presented in this paper is listed in Table 
[1] Our models are intended to be generic, but, for deflnite- 
ness, we choose the cluster A2199 as our fiducial cluster. We 





Figure 2. Time evolution of gas temperature (upper panel) and 
mass accretion rate (lower panel) at r = 5 kpc for runs A [solid 
line), Bl (short- dashed line), B2 (long-dashed line), and B3 (dot- 
ted line). See Table [T] for additional information. 



compare the steady-state profiles of electron number density 
and temperature with the observational profiles as well. 



4.1 Simulation setup 

We u se the ZEUS-3D hydrodynamic code (jStone fc NormanI 
1 19921 ) in its one-dimensional mode; we gratefully acknowl- 
edge Mateusz Ruszkowski for supplying us with the modified 
version described in RB02, which includes radiative cooling 
and thermal conduction. We solve equations ((13} , (|28|l , p9|) 
and (|30p for our fiducial cluster Abell 2199; in particular, 
we have incorporated into ZEUS a background gravitational 
potential (eq. [31]), cosmic-ray heating, cosmic-ray pressure 
support, cosmic-ray transport and the cosmic-ray energy 
equation (eq. [13]). For numerical stability, the conduction 
term is integrated using time steps that satisfy the Courant 



0.1 



a 

o 



0.01 



0.001 



"TTT 1 r 



TTT 3 




(kpc) 



Figure 3. Final density and temperature profiles after a Hubble 
time for runs A (solid line), Bl (short-dashed line), B2 (long- 
dashed line), and B3 (d otted line); crosses indicate Chandra data 
ll Johnstone et al.|[200a) . None of these models provide an ade- 
quate fit to the data. 



condition 



1 EjAr)'' 
At s; --^ '—. 

2 fKSpT 



(38) 



The time steps in our simulation are also chosen to be small 
enough to satisfy the Courant conditions required by numer- 
ical stability of the cosmic-ray energy equation 

At s; min ( — ^^=^ ) . (39) 



Our computational grid extends from rmin ~ 1 kpc to 
'"max ~ 200 kpc. In order to resolve adequately the inner 
regions, we adopt a logarithmically spaced grid in which 
(Ar)i+i/(Ar)i = (rmax/rmin)^''^, whcrc N is the number of 



active zones. We performed our main simulation (run C) in 
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three diflFerent resolutions: N = 100, 200, 400. The results of 
these three simulations are quite similar, with virtually iden- 
tical results for the second two. Therefore, we are confident 
that our simulations are numerically convergent. The stan- 
dard resolution of our simulations presented in this paper is 
iV = 400. 

For initial conditions, we assume the ICM to be isother- 
mal at T = 4.6 keV, and solve for hydrostatic equilibrium. 
We assume that at the outer boundary rmax, ?ie(?'max) = 
0.0015 cm^'^, which is close to the value extrapolated 
from the observational density profile. We assume that the 
cosmic-ray energy density Ec is a very small constant value 
throughout the cluster initially. Our results are not sensitive 
to this value, which is chosen to be — 3.8 x 10"^'* ergs s~^ 
cm""' in the models presented in the rest of the paper. For 
boundary conditions, we assume that the gas is in contact 
with a thermal bath of constant temperature and pressure 
at the outer radius, where the cooling time exceeds the Hub- 
ble time. Thus, we ensure that temperature and density of 
the thermal gas at the outer radius are constantlfl We ex- 
trapolate all hydrodynamic variables from the active zones 
to the ghost zones by allowing them to vary as a linear func- 
tion of radius at both the inner and outer boundaries. The 
intracluster gas is allowed to flow in and out of active zones 
at both the inner and outer boundaries. Cosmic ray injec- 
tion by the central AGN is only allowed when the gas at the 
inner boundary flows inward. 

4.2 The cooling flow model without any heating 

To establish a control, our first simulation (run A) follows 
the evolution of the ICM under the pure radiative cooling 
for a Hubble time t-a = {Ho = 70 km s~^ Mpc"^). The 
evolution is simple. Figure [2] shows the time evolution of the 
gas temperature and mass accretion rate at r = 5 kpc. As is 
clearly seen, the ICM cools catastrophically (characterized 
by a rapid decrease in the central temperature) and finally 
reaches a quasi-steady state, where a strong cooling flow 
(M ~ -200 Mq yr"^) is formed. Note that the long-term 
evolution of the ICM with a strong cooling flow in our model 
may not be accurate: for instance, we neglect the deepening 
of the gravitational potential well due to the large mass de- 
posit ion at the cluster center by the cooling flow fsee lMeiksinI 
Il990l ). But this control establishes a minimal baseline for the 
amount of cooling expected if there are no heating sources. 

4.3 The model with conduction only 

In this subsection, we explore the role of thermal conduc- 
tion in the evolution of the ICM. We flrst performed two 
simple simulations with / = 0.4 (run Bl) and / = 0.8 (run 
B2), following the evolution of the ICM subject to radia- 
tive cooling and thermal conduction. The time evolution of 
the gas temperature and mass accretion rate at r = 5 kpc is 
shown in Figure[2l while Figure |3]shows the final density and 
temperature profiles as a function of radius. Agreeing with 
[Cactz (1989), run Bl shows that a moderate level of ther- 
mal conduction delays the cooling catastrophe and reduces 

^ These are the same boundary conditions used by RB02. 



the mass accretion rate at the final quasi-steady state. How- 
ever, weak conduction cannot suppress the cooling catas- 
trophe sufficiently. Even with / = 0.4, the final mass ac- 
cretion rate at r = 5 kpc is around M ~ —70 Mq yr~^ 
and the final temperature at r = 5 kpc is around 1.1 keV, 
which is smaller than the observed value (~ 2 keV). Inter- 
estingly, this value of / = 0.4 is the valued needed to build 
an equilibrium model in which co nduction balances cooling 
(e.g. Zakamska fc NaravanI l|2003l )): we have verified this by 
building such a model which matches the observed tempera- 
ture and density profiles for A2199, and found the eigenvalue 
/ = 0.43. If one perturbs around the equilibrium state, the 
global thermal instability of such models are also claimed to 
be dynamica lly unimportant (the in stability growth time is 
~ 2 - 5 Gyr; iKim fc Nara"vanll2003ah . However, if one starts 
far from equilibrium, then evolution toward the equilibrium 
profile is not guaranteed. We explore this and related issues 
in a forthcoming paper (Guo et al 2007, in preparation). 
On the other hand, strong conduction successfully prevents 
the cooling catastrophe, just as shown in Figure (2] but the 
temperature does not drop significantly toward the cluster 
center, in violation of the observed temperature gradient. 

Nonetheless, it may be possible to fine-tune the value 
of / so that the final state (at t ~ tu) of the pure conduc- 
tion model produces a reasonably good fit to the observa- 
tional data. We thus performed a pure conduction simula- 
tion with / — 0.6 (run B3), which produced a somewhat 
better but still unsatisfactory fit to the data. With sufficient 
diligence it might be possible to find a satisfactory model, 
but the amount of fine-tuning seems excessive. Clearly, re- 
sults depend sensitively on the assumed value of /: if / is 
too large, the temperature profile is too close to isothermal; 
if it is too low, a strong cooling flow develops. Since nei- 
ther nearly isothermal nor strong cooling flow clusters are 
observed, if only conduction balances cooling then / must 
be restricted to a narrow range. Yet, the value of / required 
to explain observed temperature and density profiles profiles 
differs from cluster to cluster, and a physical explanation of 
how / self-adjusts in each cluster is missing; furthermore, a 
significant fraction of observed cluster s cannot be fit at all by 
conduction only models with / 5S 1 (|Zakamska fc NaravanI 
[2003iy By contrast, we find in i]4.4l that if a secondary heat- 
ing mechanism such as cosmic-ray heating is included, this 
fine-tuning problem is eliminated, and a broad range of / 
is permissible, with the remainder of the heating being sup- 
plied by CR heating in a self-regulating fashion. Note also 
from Figure [5] that the ICM in run B3 takes a long time 
(comparable to the Hubble time) to reach a steady state. 
As we shall see in the next subsection, by including a phys- 
ically motivated feedback heating term, the ICM will relax 
more quickly to a steady state, which produces a much bet- 
ter fit to the observational data as well. 

Our conduction model (eq. [36]) in this paper is some- 
what idealized. In reality, both electron conduction and tur- 
bulent mixing may contribute to hea t transport in c lusters 
(for a comprehensive discussion, see iLazarianl 1200^ . Note 
that the turbulent mixing model will probably suffer a simi- 
lar fine-tuning proble m (fine-tuning of the mix ing parameter 
may be required, see lKim fc NaravanI (|2003bl )). 
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Figure 4. Evolution of gas temperature [upper left), electron number density [upper right), cosmic-ray energy density [lower left) and 
accretion rate plotted as a function of distance from the cluster center for run C. The cluster relaxes to a quasi-equilibrium state at 
t ^ 0.4tHi except that the mass accretion rate takes slightly longer to adjust s to steady state (M a t t = OStn is plotted additionally to 
show its steady profile). Crosses in the upper panels indicate Chandra data (I Johnstone et al"]|2002l) . The short-long dashed curve in the 
lower left panel shows the final gas energy density. 



4.4 The model with conduction and cosmic-ray 
feedback heating 

The main results of this paper are presented in this sub- 
section. We consider the evolution of the ICM subject to 
radiative cooling, thermal conduction and feedback heating 
by cosmic rays. To illustrate our results, we present one rep- 
resentative model (run C) with the following parameters: 
/ = 0.3, e = 0.003, u = 0.3. We follow cluster evolution 
for a Hubble time tn- Figure |4] shows the radial profiles of 
gas temperature, electron number density, cosmic-ray energy 
density and mass accretion rate for a few time epochs. This 
model settles down to a stable steady state at t ~ OAtii- 
Figure [5] shows the gas temperature, electron number den- 
sity, cosmic-ray energy density and mass accretion rate as a 
function of time (in units of the Hubble time) for different 
distances from the cluster center r = 5, 10, 20, 40, 80, 160 
kpc. In the case of electron number density and cosmic-ray 
energy density, the above sequence of r corresponds to the 
curves from top to bottom. For the gas temperature, the 



trend is the opposite. In the case of accretion rate, the am- 
plitude of oscillations increases with r. These oscillations are 
caused by sound waves, which propagate across the cluster 
as it adjusts to changing boundary conditions (see RB02). 
The precise character of these sound waves depends on the 
resolution and boundary conditions. After the cluster relaxes 
to the quasi-steady state, similar numerical oscillations with 
a small amplitude and a short wavelength also appear near 
the cluster center, as is readily seen in the time evolution 
curves of cosmic ray energy density and mass accretion rate 
in Figure [5] as well as spatially in the pressure gradient and 
heating terms in Figure |6l 

The evolution of the cluster is very similar to that of 
RB02, since their AGN feedback heating is also triggered 
by the mass accretion. Strong X-ray emission in the cluster 
center leads to a gradual decrease in temperature and a slow 
increase in gas density, which in turn increases the cooling 
and thus increases the accretion rate. Cosmic ray injection 
is controlled by the mass accretion rate at the cluster center. 
As the CR injection rate increases, the CR heating rate also 
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increases and hence, the cluster does not cool in a runaway 
fashion. As can be clearly seen in Figure[5l the slow evolution 
of the cluster is followed by a cooling catastrophe at t ~ 
0.25tH- Unlike the standard cooling flow models where the 
gas cools to very low temperatures and a strong cooling 
flow forms, the feedback heating mechanism in our model 
suppresses the cooling catastrophe quickly after its onset. 
The gas temperature, electron number density and cosmic- 
ray energy density then cycle up and down for several times 
when the cluster adjusts its mass accretion rate in response 
to the relative importance of heating and cooling, and are 
stabilized within ~ O.ltH. The cluster thus relaxes to a stable 
steady state. The evolution of mass accretion rate is similar. 
After oscillations through positive and negative values, the 
mass accretion rate tends to a small constant negative value. 
This value is also approximately constant at all radii, as it 
should be for a steady state cluster. In the final state, the 
mass accretion rate is about M ~ —2.3 Mq yr~^, which is 
much smaller than accretion rates inferred from the standard 
cooling flow models, and con sistent with the observ ed upper 
bounds of M <, 12Mq yr~^ (| Johnstone et al.ll2002l ). 

Thus, a cooling catastrophe is averted in our model. 



Starting from a state far from equilibrium, the cluster re- 
laxes to a sustainable and stable steady quasi-equilibrium 
state. In the upper panels of Figure |4l t he observational 
temp erature and electron density profiles jjohnstone et al.] 
I2002D are also shown. As is clearly seen, the final steady state 
of our model produces a very good fit to the observational 
profiles. In the lower left panel of Figure |4j we also show the 
radial profile of thermal energy density in the final state. 
In steady state, the ratio of cosmic-ray pressure to thermal 
pressure (Pc/Pg) is always less than 0.1 and decreases away 
from the cluster center. This is well within upper bounds 
in ne arby rich clusters of Pc/Pg ^ 20% (|Ensslin et al.l 
( 19971') ■ Virgo and Perseus clusters) and Pc/Pg ^ 30% 
( Pfrommer fc EnfilinI (|2004l '). Coma cluster). 

The upper panel of Figure [6] shows the steady-state ra- 
tios of thermal and CR pressure gradients to the gravity. 
As is readily seen, the thermal pressure support dominates 
over the whole cluster, although cosmic rays provide a small 
amount of pressure support in the central regions of the clus- 
ter (~ O.lpg in the central ~ 30 kpc). In the lower panel of 
Figure (6] we show the relative importance of various heat- 
ing mechanisms in steady state. The heating due to thermal 
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Figure 6. Radial profiles of pressure support (upper panel) and 
relative importance of various heating mechanisms (lower panel) 
in the final state of our main simulation presented in ^4.41 (run 
C). Note that the short- wavelength oscillations in the curves are 
caused by sound waves due to changing boundary conditions. See 
text for additional information. 



conduction is comparable to the cosmic ray heating, which 
is dominated by wave heating through the dissipation of the 
hydromagnetic waves, while cosmic-ray colhsional heating 
is negligible. The volume-integrated cosmic-ray heating rate 
amounts to 1.1 x 10"*^ ergs s~^ in the final state, while the 
final X-ray luminosity is 2.3 x 10*^ ergs s~^. The volume- 
integrated cosmic-ray injection rate in the final steady state 
is ~ 2.4 X lO*** ergs s~^, where around half of the cosmic-ray 
energy is transported to outer regions of the cluster (> 200 
kpc) and heats the ICM in those regions. Only a small frac- 
tion (-^ 2%) of the CR energy escapes the cluster in the form 
of gamma rays and neutrinos. 

Our models and results are very robust to the level of 



thermal conduction. Runs Dl and D2 follow the evolution 
of the cluster with a lower thermal conductivity (/ — 0.1) 
and a higher thermal conductivity (/ — 0.4), respectively. 
The steady-state profiles of gas temperature and electron 
number density for both runs are very close to those for run 
C (/ = 0.3). The temperature profiles in the final state for 
both runs are shown in Figure[7K. Figures[7}3 shows the time 
evolution of gas temperature at r = 5 kpc. As is clearly seen, 
thermal conduction delays the onset of the cooling catas- 
trophe and, with a higher level of thermal conduction, the 
cluster approaches the steady state later, which agrees with 
similar results from the pure conduction models (see 14. 3p . 

We also run similar simulations with difi'erent efficiency 
of feedback (run El with e — 0.05 and run E2 with e = 
0.0003). As shown in Figure[7l3, the cooling catastrophe hap- 
pens at almost the same time (~ 0.251^)- However, with 
higher efficiency, the cooling catastrophe is more strongly 
suppressed and the final mass accretion rate is more re- 
duced. For the very low efficiency e = 0.0003, the cooling 
catastrophe is less suppressed and the gas at inner radii can 
cool to lower temperatures and higher densities in the final 
steady state, as clearly seen in Figure[7^. The final mass ac- 
cretion rate in this low-efficiency model is about M ~ — 19 
Mq yr~^, which is still much smaller than accretion rates 
in standard cooling flow models, and marginally consistent 
with the rough observ ational bound of M <J 12M0 yr~^ 
(| Johnstone et al.ll2002l ). 

Our models are also very robust to the value of v, which 
determines the spatial distribution of cosmic-ray injection 
into the ICM. Runs Fl and F2 follow the evolution of the 
cluster with a lower value (v = 0.1) and a higher value 
= 0.7) of v, respectively. The cluster relaxes to steady 
state at almost the same time as run C — 0.3); the steady- 
state profiles of gas temperature and density are also very 
similar to those for run C. Figure |8] shows the steady-state 
radial profiles of gas temperature and ratio of cosmic-ray 
heating rate (Fcr = Fwavc +rcoii) to gas cooling rate for these 
runs. As is readily seen, for higher values of u, the cosmic 
ray injection is more centrally localized, and so is the cos- 
mic ray heating. The resulting steady state mass accretion 
rate decreases with v (M ~ —6.6, —2.3 and —1.0 Mq yr~^ 
for runs Fl, C and F2, respectively). With a combination 
of cosmic ray heating and thermal conduction, our model 
produces a reasonably good fit to observation for a broad 
range of u. Note that for f < 0.1, the PdV work during the 
bubble expansion dominates over the CR injection (see Fig. 
[T]) and thus becomes the main mechanism transferring the 
AGN mechanical energy into the thermal energy of the ICM. 
We also run a simulation (run F3) with a very high value of 
u (v = 1.5), where the cosmic rays are essentially dispersed 
into the ICM at the cluster center. As shown in Figure [S] 
the cosmic rays are unable to directly heat the ICM at large 
radii and the cluster suffers a cooling catastrophe at r ~ 4 
kpc at f ~ 0.3tH. The appropriate value of u depends on the 
bubble disruption rate and is fairly uncertain (see discussion 
in t|2.3|l . However, a very high value of u seems unlikely since 
bubbles are observed to survive out to large projected radii. 

Recent studies (IChandranI l2005l : IChandran fc DennisI 
'2006*) suggest that the cosmic ray pressure gradient may 
drive convection, if the convective instability criterion 

/i dr dr 
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Figure 7. (a) Radial profile and (b) time evolution of gas temperature at r = 5 kpc for runs C {solid line), Dl [dot-dashed line), D2 
(dotted line). El (short-dashed line), and E2 (long-dashed line). See Table[T]for additional information. 




is satisfied. In cluster cores, a strong negative cosmic ray 
pressure gradient is required to drive convection, since the 
gas temperature increases away from the cluster center. 
When v is higher, the distribution of cosmic ray injection 
is more centrally peaked, and thus the resulting cosmic-ray 
pressure gradient is more negative. We checked the convec- 
tive instability criterion (eq. [40) ) for runs Fl (v — 0.1), 
C (ly = 0.3), and F2 (u = 0.7), and found that the clus- 
ter is always convectively stable during the simulations. In 
the steady state configuration of our main model (run C), 
we find that the ratio of the left-hand side to the right-hand 
side of equation (|40|) is ~ 5—10. For run F3 (u = 1.5), where 
cosmic rays are mainly dispersed into the ICM at the clus- 
ter center, the instability criterion (eq. [40]) is easily met at 
the cluster center at the very beginning of the simulation, 
suggesting that the cluster becomes convectively unstable 
long before the onset of the cooling catastrophe. Convection 
driven by the cosmic-ray pressure gradient thus provides 
an alternative means for heating the ICM and generating 
the n eeded magnetic turbulence in cluster cores (jChandranl 
I2005D . which is obviously beyond the scope of this paper. 

The CR diffusion coefficient, Kc, in galaxy clusters is 
fairly unclear. To check the dependence of our model on 
it, we performed three additional simulations with different 
radial profiles of Kc, which has the form of equation (|17p 
with Ko = 3 X 10^'^ cm^ s"^ and 3 x 10^^ cm^ s"^ for runs 
Gl and G2 respectively. For run G3, Kc is taken to be con- 
stant throughout the cluster: = 3 x 10^* cm^ s~^. The 
steady-state profiles of the contribution of cosmic ray diffu- 
sion, -Fdifi = —K,c{dEc/dr), to CR transport for these runs 
are very different, as shown in the upper panel of Figure 
[9] The middle and lower panels show the steady-state pro- 
files of CR energy density and gas temperature. These runs 
demonstrate that the evolution of the ICM is very insensi- 
tive to the radial profiles of Hc- This is due to the fact that 
the cosmic-ray diffusion time is much longer than the gas 
cooling time (see § I2.3|l . In our model, the distribution of 
the cosmic rays in the cluster is mainly determined by the 
spatial distribution of CR injection into the ICM from the 
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Figure 8. Radial profiles of gas temperature (upper panel) and 
ratio of cosmic-ray heating rate to gas cooling rate (lower panel) 
for runs C, Fl, F2, and F3. The curves for runs C, Fl and F2 
are plotted at t = OAtu, when the cluster has relaxed to steady 
state. For run F3, the curve in the lower panel is plotted at 
t = 0.25tH! which shows clearly that the cosmic-ray heating is 
centrally peaked, resulting in insufficient heating in the outer re- 
gions of the cluster and thus leading to a cooling catastrophe 
at r ~ 4 kpc as shown in the temperature profile (upper panel) 
plotted at t = O.Stn. 
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Figure 9. Radial profiles of the contribution of cosmic-ray diffu- 
sion to CR energy flux {upper panel), cosmic-ray energy density 
{middle panel), and gas temperature {lower panel) for runs C, 
Gl, G2, and G3 at t = 0.4%, when the cluster has relaxed to 
steady state. Note that the cluster relaxes to almost the same 
steady-state temperature profiles at almost the same time for 
these runs, which assume different radial profiles for Kc- See Table 
[T]for additional information. 



rising bubbles (see §[2]3]|. Note that the diffusion is also usu- 
ally subdominant with respect to the cosmic ray advection. 
The steady-state ratio of diffusion to advection for run C 
is 0.1 in the central ~ 10 kpc and becomes negligible 
(~ 0.03) in the outer regions; run G2 is an extreme case 
where the CR diffusion becomes comparable to the advec- 
tion (see Figure [9]). 



5 CONCLUSIONS 

Awareness of the significant role cosmic rays could play in 
shaping the thermal and dynamical state of gas in galaxy 
clusters has been growing in recent years. Observations of 
diffuse radio synchrotron radiation from galaxy clusters im- 
ply that strong sources of non-thermal particles are indeed 
present. At the same time, recent studies show that ac- 
tive galactic nuclei inflate buoyant bubbles containing non- 
thermal radio emitting particles and could potentially play 
a central role in suppressing the cooling flows in cool core 
clusters. Many studies have focused on the potential dynam- 
ical effects of cosmic ray pressure support, but none have 
built successful models in which cosmic-ray heating is sign- 
ficant. In this paper, we propose a new model of AGN feed- 
back heating, in which cosmic rays produced by accretion- 



triggered AGN activity heat the ICM efficiently, with only 
a small dynamical perturbation on the ICM. 

In our model, the cosmic rays are injected into the ICM 
mainly from the rising bubbles generated by central AGN 
activity, which is treated in a time-averaged sense. We as- 
sume that the cosmic rays are injected into the ICM instan- 
taneously and neglect any time delay between central AGN 
activity and the cosmic ray injection. Such time-averaging is 
justifiable because the AGN duty cycle is much shorter than 
the gas cooling time. The cosmic rays then stream along 
the magnetic field lines in the ICM. Due to the cosmic ray 
streaming instability, Alfven waves propagating nearly in 
the direction of the CR streaming are excited and scatter 
the cosmic rays in pitch angle. These waves grow exponen- 
tially until dissipated by nonlinear Landau damping, and 
thus heat the ICM efficiently. We note that the cosmic ray 
streaming may also depend on the details of the CR scatter- 
ing by the small-scale MHD turbulence in the ICM, which is 
still poorly understood (see § I2.1.1|l . Here only Alfven waves 
self-excited by the cosmic-ray streaming instability are con- 
sidered. 

We have performed a set of one-dimensional numeri- 
cal simulations of the ICM, which is subject to radiative 
cooling, thermal conduction and cosmic ray heating. If only 
thermal conduction operates, extreme fine-tuning of the con- 
duction suppression factor / is required: if / is too low, then 
a strong cooling flow develops. If / is too high, the tem- 
perature profile becomes nearly isothermal, in contrast to 
observations where the temperature invariably declines to- 
ward the cluster center. On the other hand, once cosmic 
ray heating is including, our results are very robust to the 
level of thermal conduction: the reduced cooling flow in our 
new model automatically adjusts itself to some low value 
of the mass accretion rate, which is mainly determined by 
the value of efficiency e in equation (|18|l . Furthermore, un- 
like the conduction-only case, the conduction-|-CR heating 
case rapidly equilibrates toward a steady-state solution. For 
a representative model (run C), the ICM relaxes to a stable 
quasi-equilibrium state which is a very good fit to the ob- 
served gas temperature and density profiles. The cosmic-ray 
pressure in steady state is much less than the gas pressure 
(Pc < O.lPg), while VPc ^ O.lpg in the central ~ 30 kpc 
and becomes negligible in the outer regions, all well within 
observational constraints. 

Thus, cosmic-ray heating models are a very at- 
tracti ve alternative to mec h anical heating models 
(e.g. 'Briiggen fc Kaise ij [20021: i Ruszkowski et al.1 |2004| : 
IVerna lco & Reynolds 20061: iRevnol ds ct al. 200j) in which 
the ICM is heated by the pdV work of the expanding 
bubbles, viscous dissipation of emitted sound waves or 
mixing of the hot bubble plasma with the ICM. The 
detailed microphysics of how the latter processes take place 
has not been hammered out in detail, leaving a good deal 
of uncertainty; a definitive explanation for how energy 
is transported from the observed bubbles to the ICM in 
a distributed and isotropic fashion is still outstanding. 
Which is not to say that the cosmic-ray heating model 
presented here does not suffer from similar uncertainties: 
the details of how cosmic rays leak from the bubbles, 
and/or the rate at which bubbles are disrupted, are all 
highly uncertain. Nonetheless, the results presented here 
suggest that more detailed studies of cosmic-ray heating 
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in fully 3D, cosmo l ogical si mulations (e.S-, t h e sim ulations 
of lJubelgas etal] (|2006l ): iPfrommer et alj (120061 ) where 
most of the relevant cosmic-ray physics is already included) 
are warranted. At the same time, elucidating the details 
of bubble disruption/cosmic-ray diffusion would be very 
useful in determining whether cosmic-rays or mechanical 
processses provide a more efficient means of transporting 
heat from the bubble to the ICM. 

All of these issues may assume great urgency if GLAST 
detects the 7-ray signature from the decay of neutral pi- 
ous produced when cosmic rays collide with ICM nucleons. 
For the particul a r mo del of A2199 (redshift z = 0.0309; 
Ijohnstone et al.l (|2002l )) presented here (run C), and as- 
suming that ~ 1/3 of hadronic losses go toward produc- 
ing neutral pions and hence gamma-rays, we find a steady- 
state gamma-ray fiux of ~ 9.4 x 10~^''ergs~^ cm~^, which 
varies up to a factor of ~ 2 — 3 for other runs, depending 
on the model parameters, and should be a ~ 3(t detection 
for GLAST. Since AGN activities in real clusters are likely 
episodic, we note that this value may be viewed as a time- 
averaged estimate and the real 7-ray flux may be somewhat 
different, depending on the AGN duty cycle and the cosmic- 
ray injection rate (the maximum value of 7-ray flux dur- 
ing the cluster evolution in our main simulation (run C) is 
around twice greater than that in the final steady state). 
For fainter clusters, it should be poss ible to sta ck signals to 
provide a population-averaged limit (|Ando fc Nagaii,2007l '). 
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APPENDIX A: THE TIME-DEPENDENT 
COSMIC-RAY EQUATIONS 

In this Appendix, we derive the time-dependent energy 
equation for the cosmic rays in a backgroun d magnetiz ed 
plasma (also see iMacKenzie & Voelkl (j 19821 '): lKol ll99lf )). 
When the cosmic rays are streaming along the magnetic field 
hnes at a speed faster than the local Alfven speed, they ex- 
cite hydromagnetic waves ( e.g., Alfven w a ves) by the cosmic- 
ray s treaming instability l|Lerchel 1 19671 : iKulsrud & Pearc3 
Il969l ). The cosmic rays are then scat tered in pitch an- 
gle and confined by these waves (see ISkilling|[l97ll ). for 
which we only consider forward Alfven waves propagat- 
ing nearly parallel to the background magnetic field in 
the direction of the local cosmic- ray streaming ( ba ckward 
Alfven waves are damped; see iKulsrudI l|2005l ). iLerchd 
Hill), iKulsrud & Reared il96^)). Including a net particle 
source function Q (other than compression or expansion), 
the cosmic-ray trans port equation may be written as (e.g. 
ISkilling|[l97ll . ll975al ) 



dt 



+ {u + va) ■ V/p =V • (ftpnn ■ V fp) 



+ Ip^^-{u + va) + Q, (Al) 

where fp{x,p,t) is the isotropic part (in momentum) of the 
cosmic- ray phase space distribution function, V = d/dx, 
u is the velocity of the background plasma, va is the local 
Alfven velocity, n — va/va is a unit vector along the local 
magnetic field, and Hp is the diffusion coefficient given by 



Kp(x,p) — V 



2u{x,p,ii) 



V 

T 



u{x,p,ii) 



dii, (A2) 



where v = p[l+p^/(m^c^)]^^^^/m is the cosmic-ray particle 
speed, jj, = p • n/p is the pitch-angle cosine and v{x,p,y,) 
is the frequency of pitch angle scattering of cosmic rays by 
hydromagnetic waves. The diffusion term in equation (|A1|) . 
V-(fi:pnn -V fp), shows that the diffusion follows the mag- 
netic field lines. 

From equation HA1|) . one can easily derive 



^f^ 



1 d 



^ + V-S-Q + —, — [p\u + va)-VU]=Q, (A3) 



dt 



3p^ dp 



where we have defined the particle current (e.g. IVolk|[l983l : 
IWebb & Gleesonlll979l) 

1 d -i 

S = {u + VA)fp-tipnn-V fp- — {u + VA)-^{p fp). (A4) 

op 



Integration of equation (IA3[) over all particle momenta re- 
sults in the equation for the number density of cosmic rays 
(defined as ncn = /o°° /p47rp^dp) 



dt 



Anp^Sdp - / Anp^Qdp = 0. (A5) 



To derive the energy equation for cosmic rays, we define 
three macroscopic quantities: the cosmic-ray pressure Pc, the 
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cosmic-ray energy density Ec and the cosmic-ray energy flux 
J^c as the following moments of /p and •S': 

P-^Y Jo "^""^"'^^ ^^^^ 

Ec = i-K p^Tpfpdp (A7) 
Jo 

roc 

Fc = 47r / p^TpSdp, (A8) 
Jo 

where Tp{p) is the kinetic energy of a cosmic-ray particle 
with momentum p and mass m, 



1 + 



2 \ 1/2 
p2 \ 



(A9) 



By taking the appropriate moments of equations HA1|) 
and (|A4|) and using the definitions in equations (|A6|I . HA7|) 
and (IA8I). we obtain 



-Fc = (£c + Pc)(w + wa) - nKc(n ■ VE,), (AlO) 
^ = (u + t;A)-VPc-V-Fc + Q, (All) 
where Kc(x) represents an effective diffusion coefficient 
!^P^TpKp{x,p)(n-Vfp)dp 

and Q is the net source of mean kinetic energy density of 
cosmic rays 

r oo 

Q = An p^TpQdp. (A13) 
Jo 

The fist term in the right-hand side of equation (|Alip rep- 
resents the energy-loss rate of cosmic rays due to the work 
of cosmic-ray pressure gradient on the background plasma 
(u ■ VPc) and the generati on of Alfven waves {va • V Pc, see 
iMacKenzie fc Voelkl (|l982l )). 

Assuming that cosmic rays are ultra-relativistic (i; ~ c), 
from equations HA6|I and (|A7|) . we find Pc = (7c — 1)-E'c, 
where 7^ = 4/3 is the adiabatic index for the cosmic rays 
(7c — 5/3 for non-relativistic particles). From equations 
(|A10|) and (|A11|) . we thus obtain the time-dependent cosmic- 
ray equations 

-Fc = 7cSc(u + «a) - nKc(n • V£c), (A14) 
BE 

^ = (7c - 1)(w + «a) -VBc - V ■ Fc + Q. (A15) 



